# Utilizar la librería skimage.io para leer la imagen 'histo_x.jpg' en formato RGB.
# Normalizar la imagen para que los píxeles se encuentren en el rango [0, 1]
# Visualizar la imagen
import skimage.io as io
import matplotlib.pyplot as plt

image = io.imread('imgs/histo_1.jpg')

print("Valor máximo antes de normalizar",image.max())

image = image / 255

print("Valor máximo después de normalizar",image.max())

print(image.shape)

plt.imshow(image)
plt.show()

Valor máximo antes de normalizar 255
Valor máximo después de normalizar 1.0
(1024, 1024, 3)

# Extraer la componente magenta de la imagen (que corresponde a la región tisular)
# Visualizar la imagen del canal magenta

#Al ser RGB, el magenta es la combinación de rojo, azul y un poco de verde

rojo = image[:,:,0]
verde = image[:,:,1]
azul = image[:,:,2]

mascara = (rojo > 0.5) & (verde < 0.5) & (azul > 0.5)

magenta = image.copy()

magenta[~mascara] = 0

plt.title("Magenta")
plt.imshow(magenta, cmap='gray')
plt.show()

def pink_extraction(image):

    rojo = image[:,:,0]
    verde = image[:,:,1]
    azul = image[:,:,2]

    mascara = (rojo > 0.5) & (verde < 0.5) & (azul > 0.5)

    magenta = image.copy()

    magenta[~mascara] = 0

    plt.title("Magenta")
    plt.imshow(magenta, cmap='gray')
    plt.show()
    return magenta

magenta = pink_extraction(image)

#Extracción directa de treshold
import cv2

image_uint8 = (image * 255).astype('uint8')

img = cv2.cvtColor(image_uint8, cv2.COLOR_RGB2GRAY)

plt.title("Histograma")
plt.hist(img.flatten(), bins=256, range=[0,256])
plt.show()

plt.title("Imgagen en escala de grises")
plt.imshow(img, cmap='gray')
plt.show()

mascara_magenta = (img > 200).astype('uint8')
mascara_magenta = ~mascara_magenta

plt.title("Máscara")
plt.imshow(mascara_magenta, cmap='gray')
plt.show()

# Aplicar un filtro gaussiano de tamaño 5x5 y después utilizar el método de Otsu de manera que
# los píxeles correspondientes al lumen y al background de la imagen sean 1s y el resto de los píxeles tengan un valor de 0.
# Nota: Recordar que el método de Otsu requiere como input una imagen en el rango [0-255] en formato "uint8".
# Visualizar la máscara resultante

from cv2 import GaussianBlur

def zone_extraction(img_path):
    image = io.imread(img_path)

    plt.title("Imagen original")
    plt.imshow(image)
    plt.show()

    image_uint8 = (image).astype('uint8')

    img = cv2.cvtColor(image_uint8, cv2.COLOR_RGB2GRAY)

    gauss = GaussianBlur(img, (5, 5), 0)

    plt.title('Filtro Gaussiano')
    plt.imshow(gauss, cmap='gray')
    plt.show()

    t, num_mask = cv2.threshold(gauss, 0, 1,cv2.THRESH_OTSU)
    #Se convierte a boleano
    mask = num_mask == 1

    inv_mask = (~mask * 255).astype("uint8")


    plt.title('Máscara Otsu t=' + str(t))
    plt.imshow(mask, cmap='gray')
    plt.show()

    filter_img = cv2.bitwise_and(image, image, mask=inv_mask)
    plt.title("Imagen original con filtro inverso")
    plt.imshow(filter_img)
    plt.show()

    return mask, img, image

mask, img, original_image = zone_extraction("imgs/histo_1.jpg")

# Utilizar la librería skimage.morphology.remove_small_objects para eliminar aquellos objetos cuya área sea menor a 300 píxeles
# Más información en https://scikit-image.org/docs/dev/api/skimage.morphology.html#skimage.morphology.remove_small_objects
# Visualizaer la máscara resultante

from skimage import morphology


mask_small = morphology.remove_small_objects(mask, 300)

#Aplicamos la máscara a la imagen
plt.title("Eliminación de objetos pequeños")
plt.imshow(mask_small, cmap="gray")
plt.show()

# Aplicar el algoritmo de expansión a partir de semillas (region growing) de manera que únicamente los lúmenes sean blancos
# y el resto de la imagen negra. Pista: utilizar dos semillas. Nota: Se pueden fijar las semillas de manera manual, pero
# se valorará positivamente a aquell@s que desarrollen una función para encontrarlas automáticamente.
# Visualizar la máscara resultante.
from skimage import measure
import numpy as np

def find_seeds(mask, min_area):

    labeled = measure.label(mask * 1)
    labeled = cv2.convertScaleAbs(labeled)

    seeds = []
    contours, _ = cv2.findContours(labeled, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

    #De esta forma obtenemos los centroides de las regiones identificadas.
    areas = []
    for contour in contours:
        area = cv2.contourArea(contour)
        areas.append(area)
        #Luego de examinar las areas se agrega un filtro para solo colocar
        #semillas en los lugares deseados
        if area >= min_area:
            border_pixels = contour.squeeze()
            if len(border_pixels) > 0:
                seed = tuple(border_pixels[0])
                seeds.append(seed)
                continue


    print(sorted(areas, reverse=True))
    return seeds, labeled

def apply_flood_fill(img, seeds):

    h, w = img.shape
    ref = np.zeros((h+2, w+2), np.uint8)

    img_copy = img.copy()

    img_copy = img_copy.astype("uint8")

    # seeds = [(0, 0),(1000, 1000)]

    for seed in seeds:
        cv2.floodFill(img_copy, None, seed, 0)

    # cv2.floodFill(img_copy, None, (0, 0), 255)

    img_copy = cv2.bitwise_not(img_copy)

    return img_copy

seeds, labeled = find_seeds(mask_small, 11e4)

new_mask = apply_flood_fill(mask_small * 255, seeds)
new_mask = cv2.bitwise_not(new_mask)

img_copy = img.copy()

for seed in seeds:
    cv2.circle(img_copy, tuple(seed) , 5, 255, thickness=20)

print("Semillas encontradas en", seeds)
plt.title("Semillas encontradas")
plt.imshow(img_copy, cmap="gray")
plt.show()

plt.title("Obtención de solo los lumens")
plt.imshow(new_mask, cmap="gray")
plt.show()

plt.title("Subsecciones de la mascara")
plt.imshow(labeled, cmap="gray")
plt.show()

[316521.5, 115006.5, 12884.0, 12008.0, 6076.5, 3305.0, 2161.5, 1710.0, 1263.0, 1152.0, 1000.0, 767.5, 744.0, 590.0, 383.0, 325.5, 318.0]
Semillas encontradas en [(np.int32(975), np.int32(615)), (np.int32(0), np.int32(0))]

# Rellenar los lúmenes con la función binary_fill_holes de la librería scipy.ndimage.morphology
# Visualizar la máscara resultante
import scipy as sp

mask_no_holes = sp.ndimage.binary_fill_holes(new_mask)

plt.title("Subsecciones de la mascara")
plt.imshow(mask_no_holes, cmap="gray")
plt.show()

def lumen_separation(mask, min_area, img):

    mask_small = morphology.remove_small_objects(mask, 300)

    seeds, labeled = find_seeds(mask_small, min_area)

    new_mask = apply_flood_fill(mask_small *255, seeds)
    new_mask = cv2.bitwise_not(new_mask)

    img_copy = img.copy()

    for seed in seeds:
        cv2.circle(img_copy, tuple(seed) , 5, 125, thickness=20)

    print("Semillas encontradas en", seeds)
    plt.title("Semillas encontradas")
    plt.imshow(img_copy, cmap="gray")
    plt.show()

    mask_no_holes = sp.ndimage.binary_fill_holes(new_mask)


    plt.title("Eliminación de objetos pequeños")
    plt.imshow(mask_small, cmap="gray")
    plt.show()

    plt.title("Obtención de solo los lumens")
    plt.imshow(new_mask, cmap="gray")
    plt.show()

    plt.title("Subsecciones de la mascara")
    plt.imshow(mask_no_holes, cmap="gray")
    plt.show()

    return  mask_no_holes

mask_no_holes = lumen_separation(mask, 11e4, img)

[316521.5, 115006.5, 12884.0, 12008.0, 6076.5, 3305.0, 2161.5, 1710.0, 1263.0, 1152.0, 1000.0, 767.5, 744.0, 590.0, 383.0, 325.5, 318.0]
Semillas encontradas en [(np.int32(975), np.int32(615)), (np.int32(0), np.int32(0))]

# Dibujar los contornos de los lúmenes en color verde sobre la imagen original RGB. Nota: Utilizar los flags necesarios
# para que los contornos en verde sean perfectamente visibles. 
# Visualizar la imagen superpuesta

mask_no_holes_u = (mask_no_holes * 255).astype("uint8")


canny = cv2.Canny(mask_no_holes_u, 255/3, 255)
kernel = np.ones((2,2), 'uint8')
canny = cv2.dilate(canny, kernel, iterations=3)

#Mejoramos la calidad de los bordes


contours, _ = cv2.findContours(canny, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

image_with_contours = cv2.drawContours(image_uint8.copy(), contours, -1, (0, 255, 0), 2)

plt.imshow(image_with_contours)
plt.title('Image with Contours')
plt.show()

# Determinar cuál es el lumen de mayor área y hacer un crop del mismo sobre la imagen original RGB.
# Visualizar el lumen cropeado.

#Se identifica el lumen más grande
largest_contour = max(contours, key=cv2.contourArea)

#Se realiza un crop rectangular de este para tener un contexto de la imagen

x, y, w, h = cv2.boundingRect(largest_contour)

cropped_lumen = image[y:y+h, x:x+w]
cropped_mask = mask_no_holes_u[y:y+h, x:x+w]

plt.imshow(cropped_lumen)
plt.title("Mayor lumen")
plt.show()

def big_lumen(mask_no_holes, img):
    mask_no_holes_u = (mask_no_holes * 255).astype("uint8")


    canny = cv2.Canny(mask_no_holes_u, 255/3, 255)
    kernel = np.ones((2,2), 'uint8')
    canny = cv2.dilate(canny, kernel, iterations=3)

    #Mejoramos la calidad de los bordes

    contours, _ = cv2.findContours(canny, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)


    image_with_contours = cv2.drawContours(img, contours, -1, (0, 255, 0), 2)

    plt.imshow(image_with_contours)
    plt.title('Contornos en la imagen')
    plt.show()

    #Se identifica el lumen más grande
    largest_contour = max(contours, key=cv2.contourArea)

    #Se realiza un crop rectangular de este para tener un contexto de la imagen

    x, y, w, h = cv2.boundingRect(largest_contour)

    cropped_lumen = img[y:y+h, x:x+w]
    cropped_mask = mask_no_holes_u[y:y+h, x:x+w]

    plt.imshow(cropped_lumen)
    plt.title("Mayor lumen")
    plt.show()

    return cropped_mask, cropped_lumen


cropped_mask, cropped_lumen = big_lumen(mask_no_holes, original_image)

# Calcular las siguientes características del crop del lumen de mayor área, redondeando su valor hasta el cuarto decimal.
# 1) Área
# 2) Área de la bounding box
# 3) Área convexa
# 4) Exentricidad
# 5) Diámetro equivalente
# 6) Extensión
# 7) Diámetro Feret
# 8) Longitud del eje mayor
# 9) Longitud del eje menor
# 10) Orientación
# 11) Perímetro
# 12) Solidez
# 13) Compacidad


area = cv2.contourArea(largest_contour)

# Calculate bounding box area
bounding_box_area = w * h

# Calculate convex hull and convex area
hull = cv2.convexHull(largest_contour)
convex_area = cv2.contourArea(hull)

# Calculate eccentricity, major axis length, minor axis length, orientation, and solidity using regionprops
labeled_mask = measure.label(cropped_mask)
props = measure.regionprops(labeled_mask)

# Find the properties of the largest region
largest_region = max(props, key=lambda r: r.area)

eccentricity = largest_region.eccentricity
major_axis_length = largest_region.major_axis_length
minor_axis_length = largest_region.minor_axis_length
orientation = largest_region.orientation
solidity = largest_region.solidity

# Calculate equivalent diameter
equivalent_diameter = np.sqrt(4 * area / np.pi)

# Calculate extent
extent = area / bounding_box_area

# Calculate Feret diameter (maximum distance between any two points on the contour)
feret_diameter = max([np.linalg.norm(p1 - p2) for p1 in largest_contour for p2 in largest_contour])

# Calculate perimeter
perimeter = cv2.arcLength(largest_contour, True)

# Calculate compactness
compactness = (perimeter ** 2) / (4 * np.pi * area)

# Print the results rounded to the fourth decimal place
print(f"Area: {area:.4f}")
print(f"Bounding Box Area: {bounding_box_area:.4f}")
print(f"Convex Area: {convex_area:.4f}")
print(f"Eccentricity: {eccentricity:.4f}")
print(f"Equivalent Diameter: {equivalent_diameter:.4f}")
print(f"Extent: {extent:.4f}")
print(f"Feret Diameter: {feret_diameter:.4f}")
print(f"Major Axis Length: {major_axis_length:.4f}")
print(f"Minor Axis Length: {minor_axis_length:.4f}")
print(f"Orientation: {orientation:.4f}")
print(f"Perimeter: {perimeter:.4f}")
print(f"Solidity: {solidity:.4f}")
print(f"Compactness: {compactness:.4f}")

Area: 15535.0000
Bounding Box Area: 39105.0000
Convex Area: 28602.0000
Eccentricity: 0.8445
Equivalent Diameter: 140.6406
Extent: 0.3973
Feret Diameter: 246.4650
Major Axis Length: 231.6881
Minor Axis Length: 124.0865
Orientation: 0.9970
Perimeter: 1117.6022
Solidity: 0.4872
Compactness: 6.3981

import math
import pandas as pd


def get_metric_data(cropped_mask):

    new_lab, new_num =measure.label(cropped_mask, return_num=True)


    # Crear un DataFrame para almacenar las características geométricas
    columns = ['ID', 'Area', 'BBox Area', 'Convex Area', 'Eccentricity', 'Equiv Diameter', 'Extent', 'Major Axis', 'Minor Axis', 'Orientation', 'Perimeter', 'Solidity', 'Compactness', 'Rectangularity']
    data = []

    for i in range(1, new_num+1):
        objeto = new_lab == i
        prop = measure.regionprops(objeto.astype(np.uint8))[0]

        features = {
            'ID': i,
            'Area': np.round(prop.area, 4),
            'BBox Area': np.round(prop.bbox_area, 4),
            'Convex Area': np.round(prop.convex_area, 4),
            'Eccentricity': np.round(prop.eccentricity, 4),
            'Equiv Diameter': np.round(prop.equivalent_diameter, 4),
            'Extent': np.round(prop.extent, 4),
            'Major Axis': np.round(prop.major_axis_length, 4),
            'Minor Axis': np.round(prop.minor_axis_length, 4),
            'Orientation': np.round(prop.orientation, 4),
            'Perimeter': np.round(prop.perimeter, 4),
            'Solidity': np.round(prop.solidity, 4),
            'Compactness': np.round(4 * math.pi * prop.area / prop.perimeter**2, 4),
            'Rectangularity': np.round(prop.area / prop.bbox_area, 4)
        }
        data.append(features)

    df = pd.DataFrame(data, columns=columns)
    return df

get_metric_data(cropped_mask)

def pipeline(img_path, min_area):
    mask, img, original_image = zone_extraction(img_path)
    magenta = pink_extraction(original_image)
    mask_no_holes = lumen_separation(mask, min_area, img)
    cropped_mask, cropped_lumen = big_lumen(mask_no_holes, original_image)
    df = get_metric_data(cropped_mask)

    return df, cropped_mask, cropped_lumen, img


df2, cropped_mask2, cropped_lumen2, img2 = pipeline("imgs/histo_2.jpg", 3e3)
df2

[843344.5, 2194.5, 1966.0, 974.0, 884.0, 874.5, 730.0, 684.0, 604.0, 587.5, 483.5, 395.0, 367.5, 342.5, 330.0]
Semillas encontradas en [(np.int32(573), np.int32(0))]

	ID	Area	BBox Area	Convex Area	Eccentricity	Equiv Diameter	Extent	Major Axis	Minor Axis	Orientation	Perimeter	Solidity	Compactness	Rectangularity
0	1	4675.0	7392.0	4875.0	0.8392	77.1518	0.6324	104.9016	57.0514	-0.5622	287.5635	0.9590	0.7104	0.6324
1	2	385.0	578.0	399.0	0.8947	22.1404	0.6661	34.8248	15.5540	-0.2476	88.6274	0.9649	0.6159	0.6661

Proyecto 1 computer vision

Miguel López

ID 1001014378

Convocatoria 1 - Proyecto 1¶

0) Cargar una de las imágenes histológicas¶

1) Realizar una transformación de color para convertir la imagen al espacio de color CMYK¶

2) Umbralizar la imagen para separar los píxeles del fondo de la región tisular¶

3) Limpiar la imagen eliminando los artefactos de lumen (objetos blancos pequeños que no son lúmenes)¶

4) Rellenar con 0s el fondo de la imagen para quedarnos únicamente con los lúmenes¶

5) Rellenar los objetos de los lúmenes¶

6) Detectar y dibujar los contornos de los lúmenes sobre la imagen original¶

7) Identificar y cropear el lumen más grande¶

8) Extraer 13 características geométricas que permitan caracterizar el lumen recortado¶

2. Pipeline para hist.2¶